clear; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NxN_sumstats_GZ.m
% 
% Code for "Distortions in Production Networks"
% Bigio and La'O 2020
% produces figures and tables of summary statistics
%
% inputs IOparams.mat, IOseries.mat, IOshock.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


load('IOparams.mat');
load('IOseries.mat');
load('IOshock.mat');

beg_year=1997;      % first year
end_year=2014;      % final year
N_years=18;         % Number of Years
NsecGZ=size(phi_GZ_mat,1); % Number of Industries
dYearsvec=beg_year+1:1:end_year;



%% Compute Across-time Variances and Correlations for each sector

% What is the variance of these variables over time?

% Standard Deviations
std_y=ones(NsecGZ,1); %real output
std_l=ones(NsecGZ,1); %hours
std_x=ones(NsecGZ,1); %input bundle
std_g=ones(NsecGZ,1); %sales
std_h=ones(NsecGZ,1); %expenditure-to-sales ratio
std_phiGZ=ones(NsecGZ,1); %phi GZ

%Correlations
corr_y_l=ones(NsecGZ,1); 
corr_y_x=ones(NsecGZ,1); 
corr_l_x=ones(NsecGZ,1); 

corr_y_phiGZ=ones(NsecGZ,1); 
corr_l_phiGZ=ones(NsecGZ,1); 
corr_x_phiGZ=ones(NsecGZ,1); 
corr_h_phiGZ=ones(NsecGZ,1); 


for ii=1:NsecGZ
    
    std_y(ii)=std(Dlogy_mat(ii,:));
    std_l(ii)=std(Dlogl_mat(ii,:));
    std_x(ii)=std(Dlogx_mat(ii,:));
    std_g(ii)=std(Dlogg_mat(ii,:));
    std_h(ii)=std(Dlogh_mat(ii,:));
    std_phiGZ(ii)=std(Dlogphi_GZ_mat(ii,:));
    
    corr_y_l(ii)=corr(Dlogy_mat(ii,:)',Dlogl_mat(ii,:)');
    corr_y_x(ii)=corr(Dlogy_mat(ii,:)',Dlogx_mat(ii,:)');
    corr_l_x(ii)=corr(Dlogl_mat(ii,:)',Dlogx_mat(ii,:)');
    
    corr_y_phiGZ(ii)=corr(Dlogy_mat(ii,:)',Dlogphi_GZ_mat(ii,:)');
    corr_l_phiGZ(ii)=corr(Dlogl_mat(ii,:)',Dlogphi_GZ_mat(ii,:)');
    corr_x_phiGZ(ii)=corr(Dlogx_mat(ii,:)',Dlogphi_GZ_mat(ii,:)');
    corr_h_phiGZ(ii)=corr(Dlogh_mat(ii,:)',Dlogphi_GZ_mat(ii,:)');
end


%Take Means across these variables
mstd_y=mean(std_y);
mstd_l=mean(std_l);
mstd_x=mean(std_x);
mstd_g=mean(std_g);
mstd_h=mean(std_h);
mstd_phiGZ=mean(std_phiGZ);

mcorr_yl=mean(corr_y_l);
mcorr_yx=mean(corr_y_x);
mcorr_lx=mean(corr_l_x);


mcorr_yphiGZ    =mean(corr_y_phiGZ);
mcorr_lphiGZ    =mean(corr_l_phiGZ);
mcorr_xphiGZ    =mean(corr_x_phiGZ);
mcorr_hphiGZ    =mean(corr_h_phiGZ);

medstd_phiGZ      =median(std_phiGZ);
medcorr_yphiGZ    =median(corr_y_phiGZ);
medcorr_lphiGZ    =median(corr_l_phiGZ);
medcorr_xphiGZ    =median(corr_x_phiGZ);
medcorr_hphiGZ    =median(corr_h_phiGZ);

table_corrGZ =  [mstd_phiGZ     mcorr_yphiGZ        mcorr_lphiGZ        mcorr_xphiGZ    mcorr_hphiGZ;...
                medstd_phiGZ    medcorr_yphiGZ      medcorr_lphiGZ      medcorr_xphiGZ  medcorr_hphiGZ];



% table in online appendix
disp('Table in online appendix');
disp('Sectoral distortions within-sector time-series statistics');
table_corrGZ'


%% Compute cross-sectional statistics for each variable

% Means
xsec_mDy=ones(N_years-1,1);
xsec_mDl=ones(N_years-1,1);
xsec_mDx=ones(N_years-1,1);
xsec_mDg=ones(N_years-1,1);
xsec_mDh=ones(N_years-1,1);
xsec_mDphiGZ=ones(N_years-1,1);

% Standard Deviations
xsec_std_y=ones(N_years-1,1); %real output
xsec_std_l=ones(N_years-1,1); %hours
xsec_std_x=ones(N_years-1,1); %input bundle (method 1)
xsec_std_g=ones(N_years-1,1); %sales
xsec_std_h=ones(N_years-1,1); %expenditure-to-sales ratio
xsec_std_phiGZ=ones(N_years-1,1); %phi




for jj=1:N_years-1
    xsec_mDy(jj)=mean(Dlogy_mat(:,jj));
    xsec_mDl(jj)=mean(Dlogl_mat(:,jj));
    xsec_mDx(jj)=mean(Dlogx_mat(:,jj));
    xsec_mDg(jj)=mean(Dlogg_mat(:,jj));
    xsec_mDh(jj)=mean(Dlogh_mat(:,jj));
    xsec_mDphiGZ(jj)=mean(Dlogphi_GZ_mat(:,jj));
    
    xsec_std_y(jj)=std(Dlogy_mat(:,jj));
    xsec_std_l(jj)=std(Dlogl_mat(:,jj));
    xsec_std_x(jj)=std(Dlogx_mat(:,jj));
    xsec_std_g(jj)=std(Dlogg_mat(:,jj));
    xsec_std_h(jj)=std(Dlogh_mat(:,jj));
    xsec_std_phiGZ(jj)=std(Dlogphi_GZ_mat(:,jj));
    
end


%normalize standard devations
multiple_xsec_y     =1/xsec_std_y(1);
multiple_xsec_l     =1/xsec_std_l(1);
multiple_xsec_h     =1/xsec_std_h(1);
multiple_xsec_phiGZ =1/xsec_std_phiGZ(1);

norm_xsec_std_y     =multiple_xsec_y*xsec_std_y;
norm_xsec_std_l     =multiple_xsec_l*xsec_std_l;
norm_xsec_std_h     =multiple_xsec_h*xsec_std_h;
norm_xsec_std_phiGZ =multiple_xsec_phiGZ*xsec_std_phiGZ;

%% Figures

% Figure in main text
% cross-sectional mean and std dev of phi_GZ and epsilon
figure; 
subplot(1,2,1);
plot(dYearsvec,xsec_mDh,'-.','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on; 
plot(dYearsvec,xsec_mDphiGZ,'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on; 
legend('\Deltalog\epsilon','\Deltalog\phi','Location','southwest');
title('cross-sectional mean');
subplot(1,2,2);
plot(dYearsvec,norm_xsec_std_h,'-.','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on; 
plot(dYearsvec,norm_xsec_std_phiGZ,'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on; 
title('cross-sectional standard deviation');

% Figure in online appendix
% cross-sectional mean and std dev of output, labor, and phi_GZ
figure;
subplot(1,2,1);
yyaxis left
plot(dYearsvec,xsec_mDy,'-.','LineWidth',3.5,'Color',[128/255 0/255 0/255]); axis tight; grid on; hold on; 
plot(dYearsvec,xsec_mDl,':','LineWidth',3.5,'Color',[0.5 0.5 0.5]); axis tight; grid on; hold on;
yyaxis right
plot(dYearsvec,xsec_mDphiGZ,'LineWidth',3.5,'Color',[25/255 25/255 122/255]); axis tight; grid on; hold on; 
legend('\Deltalogy','\Deltalogl','\Deltalog\phi','Location','southwest');
title('cross-sectional mean');
subplot(1,2,2);
plot(dYearsvec,norm_xsec_std_y,'-.','LineWidth',3.5,'Color',[128/255 0/255 0/255]); axis tight; grid on; hold on; 
plot(dYearsvec,norm_xsec_std_l,':','LineWidth',3.5,'Color',[0.5 0.5 0.5]); axis tight; grid on; hold on; 
plot(dYearsvec,norm_xsec_std_phiGZ,'LineWidth',3.5,'Color',[25/255 25/255 122/255]); axis tight; grid on; hold on; 
title('cross-sectional standard deviation');


disp('time series correlations');
corr_xsec_y_phiGZ   =corr(xsec_mDy,xsec_mDphiGZ)
corr_xsec_l_phiGZ   =corr(xsec_mDl,xsec_mDphiGZ)
corr_xsec_h_phiGZ   =corr(xsec_mDh,xsec_mDphiGZ)





%% Manufacturing: Durables and Non-Durables
% Compute statistics for Manufacturing industries
% Break down between Durables and Non-Durables
% Indices for Manufacturing sectors
Manufacturing=[8:26];
DurManu=[8:18];
NonDurManu=[19:26];



xsec_mDy_Manu=ones(N_years-1,1);
xsec_mDl_Manu=ones(N_years-1,1);
xsec_mDx_Manu=ones(N_years-1,1);
xsec_mDy_Dur=ones(N_years-1,1);
xsec_mDl_Dur=ones(N_years-1,1);
xsec_mDx_Dur=ones(N_years-1,1);
%Take Means across Manufacturing
for jj=1:N_years-1
    xsec_mDy_Manu(jj)=mean(Dlogy_mat(Manufacturing,jj));
    xsec_mDl_Manu(jj)=mean(Dlogl_mat(Manufacturing,jj));
    xsec_mDx_Manu(jj)=mean(Dlogx_mat(Manufacturing,jj));
    xsec_mDy_Dur(jj)=mean(Dlogy_mat(DurManu,jj));
    xsec_mDl_Dur(jj)=mean(Dlogl_mat(DurManu,jj));
    xsec_mDx_Dur(jj)=mean(Dlogx_mat(DurManu,jj));
end


mstd_y_manu=mean(std_y(Manufacturing));
mstd_l_manu=mean(std_l(Manufacturing));
mstd_x_manu=mean(std_x(Manufacturing));

mcorr_yl_manu=mean(corr_y_l(Manufacturing));
mcorr_yx_manu=mean(corr_y_x(Manufacturing));


%Take Means across Manufacturing: Durables
mstd_y_dura=mean(std_y(DurManu));
mstd_l_dura=mean(std_l(DurManu));
mstd_x_dura=mean(std_x(DurManu));

mcorr_yl_dura=mean(corr_y_l(DurManu));
mcorr_yx_dura=mean(corr_y_x(DurManu));


%Take Means across Manufacturing: Non-durables
mstd_y_nondura=mean(std_y(NonDurManu));
mstd_l_nondura=mean(std_l(NonDurManu));
mstd_x_nondura=mean(std_x(NonDurManu));

mcorr_yl_nondura=mean(corr_y_l(NonDurManu));
mcorr_yx_nondura=mean(corr_y_x(NonDurManu));


VarCovar_Manu_mat=[mstd_y mstd_y_manu mstd_y_dura mstd_y_nondura;...
    mstd_l mstd_l_manu mstd_l_dura mstd_l_nondura;...
    mstd_x mstd_x_manu mstd_x_dura mstd_x_nondura;...
    mcorr_yl mcorr_yl_manu mcorr_yl_dura mcorr_yl_nondura;...
    mcorr_yx mcorr_yx_manu mcorr_yx_dura mcorr_yx_nondura];

% table in online appendix
disp('Table in online appendix');
disp('Sectoral quantities within-sector time series statistics');
VarCovar_Manu_mat




